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Abstract 

As a result of multiple mini-jet production, initial conditions of the QCD plasma 
formed in ultrarelativistic nuclear collisions may be inhomogeneous, with large fluc- 
tuations of the local energy density (hot spots), and turbulent, with a chaotic initial 
transverse velocity field. Assuming rapid local thermalization, the evolution of such 
plasmas is computed using longitudinal boost-invariant 3+1-dimensional hydrodynam- 
ics. We compare the evolution in case that the speed of sound in the plasma is constant 
to one resulting from an equation of state involving a strong first order transition, with 
a minimum of the velocity of sound as a function of energy density. We find that 
azimuthally asymmetric fluctuations and correlations of the transverse energy flow, 
dE±/dydcf>, can develop in both cases due to the initial inhomogeneities. Hot spots 
also enhance significantly high-/cj_ direct photon yields. In the case with a phase tran- 
sition, the hadronization surface evolves into an unusual foam-like structure. Also 
in that case, we find that hadronization is considerably delayed relative to the ideal 
gas case, just as previous studies have found for homogeneous initial conditions. The 
time-delay signature of a rapid cross-over transition region in the QCD equation of 
state (as observable via meson interferometry) is thus found to be remarkably robust 
to uncertainties in the initial conditions in heavy- ion reactions. 
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DE-FG-02-93ER-40764. 
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1 Introduction 



At energies */s > 100 AGeV, copious mini-jet production in central nuclear collisions is 
expected Jl], [|, |3[] to be the main source of a very dense plasma of quarks and gluons with 
an initial (proper) energy density an order of magnitude above the deconfmement and chiral 
symmetry restoration scale [f|], e c ~ 1 GeV/fm 3 (T c ~ 160 MeV). A large number of observ- 
able consequences of the formation of this new phase of matter have been proposed based on 
a wide range of dynamical assumptions 0, and experiments are currently under construc- 
tion to search for evidence of that quark-gluon plasma (QGP) at the Relativistic Heavy-Ion 
Collider (RHIC) at Brookhaven. The observables include dilepton and hard direct photon 
yields, high-p^ jets and hadrons, strangeness and charmed hadron production, identical 
meson pair interferometry, and collective transverse expansion. Evidently, these and other 
proposed signatures depend sensitively on the assumed ensemble of initial conditions as well 
as on the transport dynamics through the hadronization point. 

In this paper we explore the dependence of several observables on the fluctuations of 
the initial conditions induced by mini-jet production. Our central dynamical assumption 
is that after a rapid thermalization time (~ 0.5 fm/c) the evolution of the plasma through 
hadronization can be approximated by non-dissipative hydrodynamics. We make this strong 
assumption to explore the maximal effects that the sought after thermodynamic properties 
of hot QCD matter may have on observables in ultrarelativistic heavy-ion collisions. Only 
hydrodynamics gives a direct link between observables and the fundamental QCD equation 
of state ||. Finite dissipative effects generally tend to reduce hydrodynamic signatures and 
require further knowledge of the transport coefficients in the QGP as well. Even in the ideal 
hydrodynamic limit, however, the observables are sensitive to the initial formation physics 
of the plasma. It is this aspect of the problem that we concentrate on in this paper. 

We show that in contrast to the conventional picture of QGP formation, the initial mini- 
jet ensemble is characterized by a wide fluctuation spectrum of the local energy density (hot 
spots) and of the collective flow field (turbulence). We also show that hydrodynamic evolu- 
tion of that inhomogeneous, turbulent ensemble can lead to novel observable consequences 
including azimuthally asymmetric transverse energy fluctuations and enhanced radiance of 
hard probes. An especially significant result is the demonstration that the time-delay sig- 
nature of the QCD phase transition, as found in previous studies 0, 0, §] starting from 
homogeneous initial conditions, survives this generalized class of more unfavourable initial 
conditions. Meson interferometry therefore remains one of the most robust and generic 
probes of QGP formation given present uncertainties of the initial conditions in heavy-ion 
collisions. 

Before elaborating on the nature of the inhomogeneous, turbulent initial conditions in- 
duced by the mini-jet production mechanism, we review first the (homogeneous) "hot-glue 
scenario" 0, [| assumed in many previous calculations of the observable consequences of QGP 
formation in A + A collisions. Mini jets are simply unresolved partons with moderate p± > 1 
GeV/c predicted from perturbative QCD. They are produced abundantly at collider energies 
because the inclusive cross section for gluon jets with moderate p±_ > po = 1 (2) GeV/c rises 
to a value <Jj e t{po) — 40 (10) mb at \/s = 200 GeV, comparable to the total inelastic cross 
section. Evidence for mini-jet production in pp and pp reactions at collider energies has been 
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inferred from the systematic rise with increasing of the yield of moderate-p^ hadrons, 
of the central rapidity density, of the enhanced tails of multiplicity fluctuations, as well as 
the flavour and multiplicity dependence of the mean transverse energy of hadrons (see Refs. 

B)- 

In the eikonal approximation to nuclear dynamics, the total number of mini-jet gluons 
produced in central A + A collisions is expected to be A 4 / 3 ~ 10 3 times larger than in pp 
collisions since each incident projectile nucleon interacts with ~ A 1 / 3 target nucleons along 
the beam axis. This simple geometric effect leads to a high rapidity density of mini-jet gluons 
with dNg/dy ~ 300 — 600 in Au + Au as shown in Fig. la. The curves are calculated via 
the HIJING model || with shadowing and jet quenching options turned off. Comparison of 
the curves for po = 1 and 2 GeV/c gives an indication of current theoretical uncertainties 
associated with the extrapolation from pp to AA reactions. The observed systematics of 
the pp data are best accounted for with p = 2 GeV/c in the HIJING model 11. That model 
combines an eikonal multiple collision formalism with the PYTHIA algorithm |10| to generate 
exclusive hard pQCD processes with p± > p and a variant of the LUND/DPM string 
phenomenology to model hadronization and account for the non-perturbative, low-p_|_ 



beam-jet fragmentation. Other parton cascade Monte Carlo models, such as developed in 
Ref . |12|] , using different structure functions and hard pQCD cutoff schemes, can account for 



the pp data using a somewhat lower p ~ 1.5 GeV/c. Theoretically [0], the scale separating 
the perturbative QCD domain from the non-perturbative (beam-jet) domain may be as low 
as po = 1 GeV/c, although no hadronization phenomenology has yet been developed with 
such a low scale that could account for the available data. Another source of moderate- 
Pj_ gluons in very heavy-ion reactions has recently been proposed based on a semi-classical 
treatment of the non-Abelian Weizsacker- Williams gluon fields [|13| . The above uncertainties 
in the initial conditions on the parton level are seen in Fig. lb to correspond to approximately 
a factor of two uncertainty of the transverse energy produced per unit rapidity in central 
Au + Au collisions at RHIC energies. 

Figure lc shows that the difference between the cases po = 1 and 2 GeV/c in the HIJING 
model is due to the production of approximately twice as many gluons in the moderate 
Pi_ < 4 GeV/c region for p = 1 GeV/c. (The -spectra extend to p± = because of initial 
and final state radiation associated with mini-jet production.) This difference is significantly 
smaller than the lowest order pQCD estimate would give because of the unitarized eikonal 
formalism used in HIJING to calculate multiple collisions and multiple jet production. For 
Po = 1 GeV/c the mini-jet cross section is comparable to the inelastic cross section. Due to 
Glauber multiple collision shadowing, the number of mini jets in that case must scale less 
rapidly than with the number of binary pp collisions. 

In Figure Id the hadronization mechanism of the mini-jet gluons via the string frag- 
mentation mechanism is found to approximately double the final hadron transverse energy 
distribution relative to Fig. lb. This is due to the pedestal or "string" effect and per- 
sists up to LHC energies in this model. The mini-jet gluons are represented as kinks in 
the beam-jet strings, and those kinks effectively boost the produced hadrons in the trans- 
verse direction. The difference between the mini-jet contribution (Fig. lb) and the final 
hadronic transverse energy distribution is due to the string model implementation of beam- 
jet fragmentation in HIJING. That component necessarily involves non-perturbative, \ow-p± 
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multi-particle production and is presently under experimental study via heavy-ion reactions 
at lower CERN/SPS energies (y/s = 20 AGeV) ||. While the string model provides an 
adequate phenomenology of beam-jet fragmentation at those energies, it is not obvious that 
it will continue to do so at RHIC and LHC This represents a significant source of theoretical 
uncertainty in calculating RHIC initial conditions. We will assume in this study that the 
extrapolation of the beam-jet physics via string phenomenology as encoded in the HIJING 
model holds up to RHIC energies. Possible sources of fluctuations of the beam-jet component 
due, for example, to "colour rope" formation have been explored in the past [14", [TJ], 



However, at collider energies, the consequences of fluctuations due to the dominant mini-jet 
component have not been considered previously to our knowledge. 

In the hot-glue scenario, the thermalization proper time is assumed to be a few times the 
mini-jet formation time r = h/po ~ 0.1 fm/c (our units are c = kg = 1). In fact, the initial 
pQCD mini-jet ^-distribution is not far from thermal as can be seen in Fig. lc, but it turns 
out that the gluon and quark multiplicities are below chemical equilibrium. Inelastic multi- 
gluon production processes are therefore essential to achieve rapid equilibration. Recent 



progress on radiative transport in QCD plasmas [17], [Tj§ [H| suggests that equilibration is 
possible on scales less than 1 fm/c. 

Taking longitudinal boost-invariant expansion into account and assuming a cylindrical 
distribution of matter, the proper energy density averaged over transverse coordinates at 



proper time r > To is given by the Bjorken formula [[20 

dE ± 1 
dy ttiR 2 



e r 



For dE±/dy ~ 1 TeV from Fig. Id, and R = 7 fm, this yields an order-of-magnitude 
estimate of e(r) ~ 65 [0.1 fm/c /r] GeV/fm 3 . If only the gluon degrees of freedom are 
equilibrated, then the temperature of the "hot glue" at the thermalization time is T(rth) — 
[e(rt/ l )/5.26] 1//4 ~ 555 (to/t^) 1 ' 4 MeV. The evolution of the temperature after r th of course 
depends on the equation of state as well as the assumed viscosity [^l[ . 



In this scenario, observables of the plasma phase such as thermal dileptons or photons 



can be computed as discussed in ||22|| . Transverse expansion can also be taken into account 
p3| as well as more realistic equations of state with a rapid cross-over transition region 
||. For transverse expansion, the temperature field, T(r, xjj, acquires a dependence on the 
transverse coordinates. In the hot-glue scenario azimuthal symmetry is naturally assumed 
for collisions at zero impact parameter. 

Another implicit assumption of the hot-glue scenario is that the initial fluid 4-velocity 
field, u^(x), vanishes in the transverse direction and that the initial flow field only reflects 
the Bjorken longitudinal expansion, 

u"(t, xj_, z) = {t/Vt 2 - z 2 , 0, 0, z/Vt 2 - z 2 ) . (2) 

Transverse expansion is allowed to develop in the course of subsequent evolution, but initially 
the plasma is assumed to be quiescent in the above sense. 

In this paper, we call into question the commonly assumed azimuthal symmetry and 
smooth radial profile of b = collisions and the above quiescent form of the initial velocity 
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field. In the next section, we show that the mini-jet mechanism does not support those 
simplifying assumptions unless the beam-jet component is much larger than estimated with 
HIJING. In Section 3, the hydrodynamic evolution of inhomogeneous, turbulent initial condi- 
tions is calculated and the novel type of azimuthally asymmetric transverse shock collectivity 
is discussed. In Section 4 the robustness of the time-delay signature associated with a phase 
transition is demonstrated. In Section 5, the enhanced radiance of direct photons from hot 
spots is estimated. A brief summary is finally presented in Section 6. 



2 The Inhomogeneous, Turbulent Glue Scenario 

Inhomogeneities arise in nuclear collisions as a result of fluctuations of the number of soft 
and hard QCD interactions per unit transverse area. Fluctuations of the soft beam-jet 
component have been considered before fl4|, [15], |l^], but at collider energies a new source of 
fluctuations that are induced by mini-jet production is expected to become dominant. Both 
types of fluctuations are strongly correlated as a function of transverse coordinates. In this 
paper, however, we consider fluctuations arising from only mini-jet production and treat the 
soft component as smooth background component of the plasma. 

Each nucleon in a central Au + Au collision suffers approximately A 1 / 3 ± A 1 / 6 ~ 6 ± 2 
inelastic collisions. Therefore, there are ~ 40 binary collisions per Oi n ~ 4 fm 2 . At RHIC 
energies, however, only a fraction aj et /ai n ~ 1/4 of those produce mini jets. The fluctuations 
of the mini-jet number density are substantial because A 1 / 3 remains relatively small even for 
the heaviest nuclei. In principle, two-nucleon correlations in the initial nuclei could reduce 
the above type of geometric fluctuations. However, the available data on high-i?x production 
in nuclear collisions at the SPS indicates sizable fluctuations, even beyond the independent 
nucleon gas approximation |24 . 



In addition to geometric sources of fluctuations, the broad transverse momentum spec- 
trum of mini-jet gluons in Fig. lc further enhances the fluctuations of the energy and mo- 
mentum deposited by mini jets per unit area. These two effects conspire to induce large 
fluctuations of the initial energy and momentum density at early times as we show below. 

The spectrum of hot spots can be computed from the HIJING event list of parton trans- 
verse and longitudinal momenta (p± a ,p Z a = P±a sinh y a ), and their longitudinal and trans- 
verse production coordinates, [z a = 0,x_|_ a ). The production points are taken from the 
initial transverse coordinates of the nucleon (diquark-quark string) in which the PYTHIA 
JTTJ| | subroutines of HIJING embed the gluons. 

To simplify the problem we assume Bjorken boundary conditions and neglect additional 
fluctuations along the longitudinal coordinate. The longitudinal velocity is thus assumed 
to vanish at z = 0. With this simplification, the transverse coordinate dependence of local 
energy and transverse momentum density in the z = plane can then be estimated from 

( MArX^O) ) = ? ( vL) P ~f F ^ ^ " X -M> 'W ■ < 3 > 

The above formula takes into account the free streaming of gluons in the transverse direction 
via x_|_ Q ,(r) = x ift + v ia r where vj_ Q = p± a /p a , as well as ideal Bjorken longitudinal expan- 
sion. The factor F(rpj_ a ) = (1 + [h/(Tp± a )] 2 )~ 1 is an estimate of the formation probability 
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T8| of partons with zero rapidity. High-p^ gluons are produced first and lower-p^ gluons 
later according to the uncertainty principle. Before r ~ h/p_i_ a , those p_i_ a -components of 
the radiation field and the evolving transient Weizsacker-Williams field of the passing nuclei 

I3| interfere strongly and cannot be treated as a kinetic gas of partons or a plasma. The 
exact form of F is not important except for extremely small times r < 0.5 fm/c prior to 
thermalization. 

The above expression, when averaged over transverse coordinates, reduces to the original 
Bjorken estimate, eq. (fl|) with (£(t)) — e(r). In addition, the averaged transverse momentum 
density and hence the flow field vanishes, (Mj_(r)) = 0, up to finite multiplicity fluctuations. 

To study the transverse coordinate dependence of the initial conditions given discreet 
mini-jet phase space coordinates, we must specify the transverse, Ar_i_, and longitudinal, Ay, 
resolution scales corresponding to an elementary fluid cell. The densities, coarse-grained on 
that resolution scale, are obtained from (|^) by the simple substitution 

5 (2) (x ± - x ±a (r)) S(y a ) - e ( A ^/ 2 -M) jj (Ar ± /2 - \x ia (r) - x<\) . (4) 

Ar I Ay 

_L a i=x,y 

The uncertainty principle limits how small Ar± > h/Ap± could be. However, in our 
dynamic problem, causality restricts how large it can be. At a given proper time r, after 
the two highly contracted nuclei pass through each other, the local horizon for any gluon in 
the comoving (y a = 0) frame only has a radius cr. Thus, at the thermalization time, r t h, 
each gluon can be influenced by only a small neighbourhood of radius cr t h — 0.5 fm of the 
plasma. In order to minimize fluctuations, we will take the transverse resolution scale to be 
the maximal causally connected diameter, Ar± = 2cr t h — 1 fm. Also, we take the rapidity 
width to be Ay = 1, since gluons with larger relative rapidities tend to be produced too 
late (due to time dilation) to equilibrate with gluons with y = 0. Note that the number, 
N(t) = [R/(ct)} 2 , of causally disconnected domains in a nuclear area, nR 2 , is initially very 
large and even at the onset of hadronization, ~ 3 fm/c, several disconnected domains 
remain. 

In Figure Id we saw that at y/s = 200 AGeV, the HIJING model predicts that approxi- 
mately half of the final produced transverse energy arises from beam-jet string fragmentation. 
As emphasized before, it is unknown whether the string model coupling between the beam- 
jet fragmentation and mini-jet fragmentation, that accounts so well for pp multi-particle 
data IJ, can be applied to collisions of heavy nuclei. Effective string degrees of freedom may 



be quenched due to the colour conductivity of the plasma |25fl , but soft field fluctuations 
associated with beam-jet fragmentation are likely to persist and contribute a background 
source of \ow-p± quanta at some level. In the present study, we do not hadronize the par- 
tons via this default JETSET string mechanism in HIJING but only use their phase space 
density to compute the mini-jet contribution to the initial energy-momentum tensor T IJ,u {x) 
as above. Since our dynamical assumption is the applicability of hydrodynamics, we model 
the soft beam-jet component by a homogeneous background fluid. We estimate the energy 
density of that component from the HIJING final hadronic dE±/dy. One advantage of the 
hydrodynamic approach is that by treating the beam-jet component as a fluid, we need not 
specify further details of its uncertain microscopic nature. In our case, the main function 
of that background fluid is to reduce the energy density fluctuations induced by the mini 
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jets. The transverse coordinate distribution of that soft component is assumed to reflect the 
density of wounded nucleons, i.e., it is taken proportional to (1 — [r±/ R] 2 ) 1 / 2 . 

If we take into account the relatively long formation time of soft, (p±) ~ 0.3 GeV/c, 
partons, only about half of the dE s °^ 1 /dy ~ 0.5 TeV from Fig. Id contributes to the soft 
background. In that case, this soft component adds a relatively small, smooth contribution 
~ 5 GeV/fm 3 to the central energy density at r th = 0.5 fm/c. 

In Figures 2a,b the initial energy and momentum density profile of a typical central 
Au + Au event is shown as a function of the transverse coordinates in the z = plane. 
This corresponds to a HIJING event with the default mini-jet cutoff scale po = 2 GeV/c. 
The profile for po = 1 GeV/c looks similar except the energy density scale increases by 
approximately a factor of two. The striking feature illustrated by the lego plot in Fig. 2a is 
the existence of several prominent "hot spots" with £ > 20 GeV/fm 3 separated by ~ 4 — 5 
fm. In this event the hottest spot reaches an energy density of about 40 GeV/fm 3 . Between 
the hot spots are cooler regions with energy density down to the soft background scale of 
~ 5 GeV/fm 3 . 

The turbulent nature of the initial conditions is illustrated in Fig. 2b. An arrow plot 
representation of the transverse momentum field is shown. The highest momentum regions 
tend to coincide with the regions of highest energy density. The initial transverse flow 
velocities are found to be distributed broadly up to ~ 0.5 c. We note that turbulence here 
is not inherent to QGP evolution since the Reynolds number of the QGP, IZe ~ R/l ~ 10, 
is not large. Thus, laminar flow is not expected to break up into complex vortex structures. 
In our case, the turbulence of the QGP initial conditions is induced by the external mini-jet 
production mechanism. This type of turbulence is analogous to flow induced by the blades 
of a mixer in a bowl of liquid. 

In Figure 2c the distribution of energy densities is shown, coarse-grained over 1 fm 3 cells 
and averaged over 200 HIJING events. Only cells in the central region with r±_ < 4 fm are 
considered in this histogram to reduce fluctuations from the less interesting surface region. 
The event-averaged energy density (£(r th )) ~ 12 GeV/fm 3 includes the soft background 
contribution discussed above. However, the distribution is highly asymmetric with relatively 
high probability of large and small fluctuations. The rms relative fluctuation of the initial 
energy density is found to be AS/ (£) ~ 0.7 for this assumed soft background level. 

In Figure 2d the distribution of the effective local gluon temperature, T eff (r t / l , xj_) = 
(S(r t h, xjJ/5.26) 1 / 4 , corresponding to Fig. 2c is shown. (This estimate of the temperature 
neglects collective flow velocities and that the gluon number as computed in HIJING is 
not in chemical equilibrium.) The local temperature is seen to fluctuate around the mean 
~ 350 MeV with an rms width of AT ~ 60 MeV. It is clear from the above that the mean 
values of £ and T do not provide an adequate characterization of such an ensemble of initial 
conditions. The ensemble of mini-jet initial conditions has a broad fluctuation spectrum 
primarily because the local event horizon is so small at early times. Fluctuations of £ in finite 
volumes arise even in an ideal Stefan-Boltzmann thermal ensemble with (e) = KT A . A simple 
estimate of the magnitude of thermal fluctuations is given by Ae/(e) ~ 2/(KT 3 V) 1 / 2 ~ 0.5 
for T = 350 MeV, K ~ 5.26 (if only the gluon degrees of freedom are assumed to be 
equilibrated), V = A(cr t h) 3 (for Ay = 1). This is comparable to the fluctuations induced 
dynamically by the mini-jet formation. The spectrum of fluctuations differs from an ideal 
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thermal one because they are driven here by Glauber nuclear geometry and the pQCD 
mini-jet spectrum. 

In Figure 3a the smooth, event-averaged energy density profile is compared to the fluc- 
tuating energy density profiles of three other separate HIJING events. The azimuthally 
symmetric, event-averaged surface corresponds to the one usually assumed in the hot-glue 
scenario. The three individual events in parts b-d, on the other hand, show that the dynam- 
ical fluctuations extend up to 40 GeV/fm 3 (see Fig. 2c) and cause strong deviations from the 
smooth, event-averaged profile. The shaded contour plots above the surface plots provide 
another representation of the azimuthally asymmetric inhomogeneities caused by mini-jet 
production at the event-by-event level. 

It is important to note that the hot spots are not due to isolated hard pQCD jets, but 
rather to the accidental coincidence of several moderate-p_i_ (see Fig. lc) mini jets in the 
same 1 fm 3 volume. This is seen by comparing the gluon number density profiles in Fig. 4 
to the corresponding energy density profiles in Fig. 3. The typical gluon density is seen to 
be about 5-10 fm~ 3 at this early time r t h = 0.5 fm/c. This includes, in addition to mini 
jets, softer gluons from initial and final state radiation as well as the soft partons from the 
beam-jet fragmentation component. 

Again we emphasize that Figs. 2-4 correspond to the minimal fluctuations from the 
point of view of hydrodynamic evolution. Only at later times, when the energy density 
is significantly lower, can the fluctuations be significantly reduced by coarse graining on 
much larger resolution scales. The dilemma is that if indeed the local thermalization is 
so short as hoped for in the hot-glue scenario, then inhomogeneous and turbulent initial 
conditions must be considered when calculating signatures of QGP formation. On the other 
hand, if the thermalization time turns out to be much longer than current estimates, then 
hydrodynamics is of course inapplicable and the observables sensitive to the initial "hot" 
phase of the reaction will not provide information about the sought after thermal properties 
of the dense QGP. For our purpose of exploring optimal signatures associated with the 
thermal properties of dense matter, we must rely on rapid thermalization and therefore live 
with inhomogeneities in calculating plasma observables. Inevitably this means that many 
observables will depend strongly on the precise form of the ensemble of initial conditions. The 
mini-jet initial conditions with a string model beam-jet background represent a particular 
ensemble obtained in extrapolating present pp phenomenology to nuclear collisions via the 
HIJING model. We will refer to the above ensemble of initial conditions as the "turbulent- 
glue scenario" to contrast it with the simpler initial conditions assumed in the hot-glue 
scenario. 

Examples of plasma probes sensitive to the highest temperature phase are high-mass 
dilepton pairs, direct photons, and heavy quark production. They are exponentially sensi- 
tive to the fluctuation spectrum of local temperatures because of the Boltzmann suppres- 
sion factor, exp(— M±/T). The hot-glue and turbulent-glue scenarios differ considerably in 
the predicted yields of such observables because the ensemble average of hydrodynamically 
evolved turbulent initial conditions and the hydrodynamic evolution of ensemble-averaged 
initial conditions are not equivalent. 
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3 Hydrodynamic Evolution of the QGP 



3.1 The SHASTA Algorithm and Tests 

In order to explore how such turbulent plasmas may evolve, we solve numerically the hydro- 
dynamic equations given the above ensemble of initial conditions. Hydrodynamic evolution 
of homogeneous plasmas has already been studied extensively || p3fl . In most previous 
studies, azimuthal (cylindrical) symmetry was assumed to simplify the problem. One excep- 
tion is the study of homogeneous but azimuthally asymmetric initial conditions produced in 
non-central nuclear collisions PE| . Our study focuses exclusively on b = collisions, where 
cylindrical symmetry of ensemble-averaged observables must always hold. At the event-by- 
event level, azimuthal symmetry is of course broken by the probabilistic nature of hadronic 
interactions. However, we are not interested in statistically uncorrelated fluctuations, rather 
only in the azimuthally asymmetric correlations that can evolve dynamically out of specific 
inhomogeneous mini- jet initial conditions. 

We assume longitudinal boost-invariance, i.e., inhomogeneities in the rapidity dimension 
are neglected in this first study. Mini jets from HIJING actually lead to hot-spot fluctuations 
that extend only one unit in rapidity. The true 3-dimensional inhomogeneities are therefore 
significantly larger than what will be studied here assuming that the fluid density is constant 
along fixed proper time surfaces r = \ft 2 — z 2 . The full treatment of the problem will require 
much more elaborate 3-dimensional simulations in the future. At present, a 3-dimensional 
code is in the development phase. 

A very important advantage of using a hydrodynamic approach is that hadronization 
can be taken into account using the non-perturbative equation of state, p(e), as deduced 
from lattice QCD simulations. The price paid for hydrodynamics is of course the necessary 
assumption that the equilibration rate is large compared to space-time gradients of the fluid 
T^-field. 

In order to compute the evolution of turbulent initial conditions we solve the equations 
of relativistic hydrodynamics, 

d^ = 0. (5) 

We rely on the most recent optimistic estimates [[17], HH, ^ that suggest that the local 
equilibration time in a QGP may be short enough to neglect dissipative effects. In that case, 
T^ u = (e + p)u^u u — pg^ u is the energy-momentum tensor for an ideal fluid, where e is the 
(proper) energy density and p the pressure in the local rest frame of a fluid element moving 
with 4-velocity = 7(1, v) in the computational frame (v is the 3-velocity, 7 = (1— v 2 ) -1 / 2 , 
and g^ u = diag(+, — , — , — ) is the metric tensor). The equations ([5]) are closed by specifying 
an equation of state p(e). 

Since we assume boost-invariance in z— direction, it suffices to solve the equations of 
motion at z — 0. Furthermore, since boost-invariance implies v z = z/t, the four eqs. (|5|) can 
be simplified to yield the three equations (valid at z = 0) 

d t £ + d x [(£+p)v*]+d y [(£+p)v y ] = -F(£,p,t), 

d t M x + d x (M x v x + p) + d y (M x v y ) = -G(M x ,t), (6) 

d t M y + d x (M y v x ) + d y (MV + p) = - G(M y , t) , 
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where F(S,p,t) = (£+p)/t, and G(M, t) = M/t. (Our notation is T 00 = £, T 0i = M\ i = 
x,y.) These equations are solved numerically via a two-step operator splitting method. 
The first operator splitting step decomposes the problem into solving the above system of 
equations with F = G = (this corresponds to purely two-dimensional fluid motion), and 
then updating the obtained solution S, M x , M y , p, ... according to the ordinary differential 
equations 

^ = -F(E,p,t), ^ = -G(M\t), i = x,y, (7) 
i.e., more explicitly one corrects 

£ = S-F(S,p,t)dt , M l = W -G(M\t)dt , i = x,y. (8) 

This method was originally suggested by Sod and was proven to be adequate for treating 
the analogous, azimuthally symmetric boost-invariant problem in [|J. 

The solution to © with F = G = itself involves another operator splitting step, i.e., 
one first solves the system of equations 

d t £ + d x [(S + p)v x ] = 0, 
d t M x + d x (M x v x +p) = 0, (9) 
d t M y + d x (M y v x ) = , 

corresponding to transport of the hydrodynamic fields in x— direction only, and with the 
solution to this set of equations one solves 

d t £ + d y [(£+p)v y ] = 0, 
d t M x + d y (M x v y ) = 0, (10) 
d t M y + d y (M y v y +p) = , 

corresponding to transport in y— direction only. Equations (|9|) and ( |10D are solved with the 
phoenical SHASTA algorithm [28] as presented in [29j], with half-step updating of the source 



terms and simplified source term treatment. The transport steps (|9|) and ( [TDD are alternated 
between successive time steps to minimize systematical errors in the propagation. Per each 
time step, the fields are propagated according to @ and fllCf ), and finally corrected with (§). 

The fluid evolution is studied for two idealized forms of the QGP equation of state. One 
is the ideal relativistic g 

p(e) = e/3 = KT 4 /3 . (11) 

Here, K ~ 5.26, corresponding to an equilibrated gluon gas. In the second case, a Bag model 
equation of state with a strong first order transition at a critical temperature T c is assumed: 

if e > €q , 

if e Q > e > e H , (12) 
if e H > e . 
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With the Bag constant B and ratio of effective plasma and hadronic degrees of freedom 
r = Kq / Kh fixed, the energy density of the mixed phase is bounded by 

4r - 1 
r — 1 

£h = — — r B . (13) 
r — 1 

We take here B = 0.39691 GeV/fm 3 and r = 37/3. For this choice, e Q ~ 1.7 GeV/fm 3 and 
e# ~ 0.1 GeV/fm 3 , the critical pressure is p c = e#/3 and the critical temperature is T c ~ 169 
MeV. Hydrodynamics with a more realistic equation of state with a finite cross-over region 
AT/T C ~ 0.1 was considered in Ref. M. In the present exploratory study, only the two 
idealized forms above will be contrasted. 

The evolution equations are solved on a two-dimensional cartesian 200 x 200 mesh with 
grid spacing Ax = 0.2 fm. The Courant-Friedrichs-Lewy number is taken as A = At/ Ax = 
0.4. The cartesian grid breaks the isotropy of space and might lead to instabilities. As 
a first check to determine whether our multi-dimensional algorithm tends to produce such 
numerical artifacts, we consider the expansion of a cylindrically symmetric Gaussian hot 
spot with radius 2 fm and peak energy density S = 30 GeV/fm 3 at rest. The time evolution 
should respect the initial cylindrical symmetry of the problem. In Figure 5a we show the 
initial (calculational frame) energy density profile (T 00 = £). Figure 5b shows the energy 
density profile after evolving the hydrodynamical equations with the standard version of 
the SHASTA algorithm as described in |2{J and the ideal gas equation of state ( |Tl| ) at 



time t = 14.9 fm/c (for the sake of clarity we show the profile for positive x, y only, the 
other quadrants are symmetric). The observed strong fluctuations which break cylindrical 
symmetry are due to the following: the flux limiter in our version of the SHASTA (cf. 
eq. (18) of f29[) prevents the occurrence of unphysical extrema only along the direction of 



propagation, i.e., the x— and y— direction. Off the grid axis, small perturbations are not 
smoothed out and can grow to produce the features seen in Fig. 5b. 

One way to cure this is to use a multi-dimensional flux limiter as proposed by Zalesak 
||30|| . Here, however, we choose the simpler method of reducing the coefficient in front 
of the antidiffusion fluxes (default value 1/8, cf. eq. (17) of El; the necessity to reduce 



this coefficient occurs also in other, purely one-dimensional situations, cf. the discussion in 
|29| , PH]). Figure 5d shows the same situation after evolving with the antidiffusion coefficient 
reduced to 70% of its default value, i.e., to 0.7/8=0.0875. This obviously strongly reduces 
symmetry breaking. Although this prescription increases the numerical diffusion of the 
algorithm and thus produces physical entropy, all our following results are generated with 
this reduced antidiffusion in order to avoid spurious azimuthal asymmetries of the transverse 
energy, dE^/dyd(j), of purely numerical origin. In Figure 5c we show the energy density after 
evolving with the equation of state with the first order phase transition ([12]) . Also in this 
case, the reduction of the antidiffusion helps to preserve the cylindrical symmetry of the 
problem. 

Another test of the algorithm is to check how well it reproduces test cases with an- 
alytically known solutions. One such test would be the expansion of a cylinder with a 
sharp surface. Here, however, we focus instead on the more physically relevant problem 
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of the expansion of the cylindrically symmetric (smooth) Gaussian studied in the previous 
figure. Although that problem is effectively one-dimensional, it does not have a purely 
analytical solution. We may, however, compare the numerical solution generated with our 
multi-dimensional SHASTA algorithm to that obtained with a one-dimensional algorithm 
which is known to reproduce analytical solutions very well, namely the relativistic Harten- 
Lax-van Leer-Einfeldt (RHLLE) algorithm [ |2"9"| , pi] , |52"| , for our purposes modified with 



a Sod operator splitting step as described in || to account for longitudinal boost-invariance 
(that operator splitting is in fact analogous to the one in eq. (B)). 

In Figure 6 we compare energy density profiles along the x— axis for (a) the ideal gas 
equation of state and (b) the Bag model equation of state with phase transition. For the 
RHLLE run we used a 2000 cell grid with Ax = 0.02/A fm, A = 0.99. The larger prediffusion 
|29| for the SHASTA visible in (a) is expected on account of the smaller Courant-Friedrichs- 



Lewy number 0.4 as compared to 0.99 used in the RHLLE run [25|. The sightly slower cooling 
of the center is due to the larger numerical diffusion of the SHASTA. The sharp cusp- 
like structures that develop in the RHLLE solution in (b) and which are associated with 



deflagration discontinuities in the transition from quark-gluon to hadron matter |Hj |4 
are therefore also broadened in the SHASTA calculation. (The numerical diffusion for the 
RHLLE is so small that it even tends to produce small-scale oscillations around the true 
solutions at the origin and for late times in (b).) Up to these small numerical effects, however, 
agreement is satisfactory and establishes confidence in that the SHASTA algorithm is able to 
generate approximately correct hydrodynamical solutions also for more complicated initial 
conditions. 

A third test of the numerical stability of the algorithm is shown in Fig. 7, comparing the 
experimentally observable transverse energy flow, dE^/dyd<f), and the azimuthal correlations 
of the transverse energy flow, Cet(A(I)), m the case °f ideal gas and Bag model equations of 
state. An initially azimuthally symmetric Gaussian energy density profile of radius 1 fm is 
evolved 100 time steps (with time step width At = X Ax = 0.08 fm/c). The thermal smeared 
transverse energy distribution is computed as a function of time via eq. ( |2"9"| ) derived in the 
Appendix. The azimuthal transverse energy correlation is defined by 

r (AM E ± (<f> + A<f>) E ± (<j>) 

c ET (^) = j o ^ 1, (14) 

where we abbreviated E±(</>) = dE±/dyd<p and (E±) = j d(p E±(4>)/2tt. For this test case, 
the exact solution has E±(</>) = (E±) independent of and Get = 0. As shown in Fig. 7, 
while the initial condition agrees with that expectation, numerical errors develop at later 
times due to the cartesian grid that breaks the symmetry of the initial conditions. This is 
especially obvious from the fact that the azimuthal anisotropy develops directed along the 
four diagonal directions of the grid. The numerical anisotropy is largest for the first order 
transition case and reaches 10% in the E±((f>) distribution while the numerical E± correlations 
remain below 0.1%. These results indicate the magnitude of numerical errors we must keep 
in mind in order to assess the significance of results for more complex geometries. 
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3.2 Transverse Shocks in Inhomogeneous Geometries 

In this section we demonstrate that a novel class of azimuthally asymmetric collective flow 
patterns can develop from inhomogeneous initial conditions. These patterns are analogous 
to the hadronic volcanoes proposed by T.D. Lee However, instead of emerging from 

instabilities at the hadronization surface, these "volcanoes" arise from transverse directed 
shock waves formed during the expansion of initial state inhomogeneities p6| . 

To illustrate this type of collectivity, consider the evolution of two idealized Gaussian hot 
spots with radii 1 fm and with centroids separated by 4 fm in the transverse x— direction. 
The initial transverse flow velocity is assumed to vanish and only the Bjorken longitudinal 
flow velocity (||) is taken into account. The initial energy density of each Gaussian is assumed 
to peak at 30 GeV/fm 3 to simulate the hot spots seen in Figs. 2-4. 

In Figure 8 the evolution of the expanding hot spots is shown for an ideal gas equation 
of state. The shaded energy density contours are shown every 10 computational time steps, 
corresponding to 0.8 fm/c between each frame. The expansion proceeds, as expected [[37] 



with the formation of two (cylindrical) shells. At 2.9 fm/c (left column, middle row), a high 
density wall of shocked matter has formed at x = where the expanding shells intersect. 
Subsequent evolution of that dense region leads to the "eruption" of two back-to-back jets 
of matter in the transverse y— direction. These transverse shock patterns are similar to the 
familiar "squeeze-out" flow pattern [38| produced in lower energy nuclear collisions. However, 



they are produced here by the collision of radially expanding shells or bubbles of relativistic 
hot matter, and in the case of multiple initial inhomogeneities, multiple volcanoes form at 
azimuthal angles that depend on the particular geometry. 

These transverse shocks are most clearly visible in the dEj_/dyd<f) distribution as shown 
in Fig. 9. The initial d£ , j_(r t / t )/d?/d0 ~ 24 GeV is of course rotation invariant. However, 
at the end of the evolution (r = 7.7 fm/c), two narrow towers of directed transverse energy 
emerge at = 90 and 270 degrees. This occurs because the initial azimuthal asymmetry 
in coordinate space is transferred into azimuthal asymmetry in momentum space through 
evolution. Finite dissipative effects would of course decrease the intensity and broaden the 
azimuthal distribution of these volcanoes. 

Note also that the overall magnitude of the transverse energy decreases with time because 
work is performed by the fluid associated with longitudinal Bjorken expansion. The Gaussian 
hot spots expand not only in the transverse direction but also along the beam direction. 

The initial ^-correlation function is of course zero. By the time of freeze-out {£ ~ 
0.1 GeV/fm 3 ), however, a strong forward and backward correlation develops for this ge- 
ometry. It is important to observe that the magnitude of the correlation is only 30% even 
though the collective signal-to-thermal noise ratio, s, in dEj_(Tf = 7.7 fm/c) /dyd(j) is peaked 
at s ~ 4. This is a general feature of correlation functions since the convolution introduces 
a dependence on the width of the signal as well. If the relative width of the signal to noise 
is S, then a simple estimate of the auto-correlation is 

C(0) - (g - 1)2 ^-^ (15) 
G(Uj " (1 + 6(8-1))* ' [ } 

We can understand qualitatively the magnitude of the correlation at A<p = in Fig. 9 taking 
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s ~ 3 and 5 ~ 0.2. The off-peak anti-correlation can also be estimated as 

( s - 

which is also in qualitative accord with the computed correlation. These correlations are 
numerically significant because from Fig. 7 we found that the numerical errors lead only to 
a 0.01% correlation in the ideal g clS CclSC. The slight forward-backward asymmetry in the 
correlation function is only due to our histogram binning procedure. 

In Figure 10 we show the evolution of the same initial condition assuming the Bag model 
equation of state, eq. ([12]). The expansion in this case is qualitatively different from the ideal 
gas case. The expanding shells or bubble walls are much thinner as is also the high density 
intersection region. This is even more clearly seen in Fig. 11 that shows the energy density 
as a function of time along the x— axis for both equations of state. The transverse velocity 
profile along that slice is also shown. In the case of a first order transition, a sharp cusp is 
produced at x ~ —5, 0, 5 fm, which correspond to points where matter cools to the critical 
mixed phase transition point, eg. In contrast, the bubbles and transverse shock zones in the 
ideal gas case remain comparable to the width of the initial hot spots. Those structures are 
much thinner in the phase transition case. The reason is that mixed phase matter does not 
resist compression due to the vanishing of the speed of sound and therefore fills less space 
than an ideal gas with finite compressibility (see also Fig. 6). 

Also, the evolution in the case of a first order transition is slower by about a factor of two 
relative to Fig. 8. This is due to the stall of the expansion in the mixed phase because of the 
vanishing velocity of sound in those regions || [33|. In the ideal gas case, we have (? s = 1/3 



throughout the expansion. With c s = in the region with e# < e < eg, sharp cusps are 
formed that move outwards only slowly. As can be seen from the velocity profiles the flow 



velocity has a discontinuity through the cusp, typical of deflagration phenomena || [33], HJ. 

The hydrodynamic stability analysis of bubble formation is complex and in the cosmo- 
logical electroweak context is still subject to controversy |39|, |40|. In the QCD case even 



less is known, though at least in |39]] marginal stability was found to be possible within the 
uncertainties of the relevant scale. In our case, bubble formation does not result from super- 
cooling but rather from the dynamical expansion of initial state inhomogeneities. Whether 
these hydrodynamic structures are stable is left here as an open question, especially since 
the thickness of the bubble walls is of hadronic dimensions. As we show in the next section 
the stalled expansion and the formation of thin shells of expanding mixed phase matter is 
the typical pattern we find also for the more complex inhomogeneous, turbulent mini-jet 
initial conditions. 

Returning to Fig. 9, we see that the consequence of stalled expansion is a considerable 
reduction of the transverse shock intensity as measured by dEj_/dyd(p, relative to the ideal 
gas case. The signal to noise is reduced to s ~ 1.5 and the relative width is also reduced 
to 5 ~ 0.1. With these reductions, the E± correlation is seen to be reduced by an order 



of magnitude in accord with eqs. (P"5"|,IT)). However, the few-percent correlation is still 
numerically significant compared to the 0.1% numerical errors deduced from Fig. 7. 
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3.3 Evolution of Turbulent Mini-Jet Initial Conditions 



In Figures 12, 13 we compare the evolution of a typical mini-jet initial condition. The initial 
rptiv required for the hydrodynamic evolution is taken from eq. (§) with T 00 = £, T 0t = M % . 
We note that the HIJING initial conditions are not of the ideal fluid form. By only taking 
the left hand side of (|3|) to fix all other components of the fluid energy-momentum tensor, 
we convert the HIJING initial conditions into an ideal fluid form through the assumption of 
thermalization at r t h- Thermalization has the effect of reducing the initial transverse energy 
somewhat from the HIJING input, because some of the transverse energy is converted into 
longitudinal thermal motion. 

The time steps between frames in the case of a first order transition (Fig. 13) are taken 
to be twice as long (1.6 fm/c) as in the ideal g (Fig. 12) to take into account the 

inherently slower and stalled expansion in the first order transition case. For this event, 
several prominent hot spot regions are seen to expand in a manner analogous to the previous 
Gaussian examples. However, in this case the hot spots also start with initial transverse 
collective flow velocities determined from the mini-jet fluctuations. 

The background fluid produced via soft beam-jet fragmentation in this example is as- 
sumed to have the smooth, azimuthally symmetric profile 

A psoft i q 

£ soft (n,,x ± ) = ^ -^^-rl/R 2 ) 112 , (17) 

dy r th TcR 2 2 

with = as well. Unlike in Figs. 2-4 where we included only about one half of the 

soft transverse energy on account of the formation time estimates, in this case we have taken 
the full dE s f ft /dy = 0.5 TeV estimated from HIJING. The soft component adds in this case 
(i.e. po = 2 GeV/c) a larger smooth background of depth £ soit (T t h, 0) ~ 10 GeV/fm 3 to the 
central region. As emphasized before, the magnitude of that soft background component is 
uncertain to about a factor two, and we take the above value to be on the conservative side. 

In Figure 12 one can see two main expanding bubbles emerging from the two hottest spots. 
Their evolution is more irregular than in Fig. 8 because of the inhomogeneous background in 
which they propagate. The solid curve shows the freeze-out hadronization contour with e = 
e H = 0.1 GeV/fni . In this ideal this hadronization surface shrinks monotonically 

but does not correspond to any obvious £ contour because of the underlying chaotic collective 
flow velocity field. 

In Figure 13, in addition to slowing down of the expansion, the evolution with the first 
order transition leads to the production of much thinner and irregular bubble fragments of 
mixed phase matter reminiscent of the thin shells in Fig. 10. It is the multiple inhomo- 
geneities and shear velocity fields that make these structures much more irregular in this 
case. The hadronization surface (e(r, x_|_) = £h) also seems to acquire a foam-like structure^ 
The hadronization surface tends in this case of a first order transition to coincide closer 
with the plotted £(t, x_|_) = T 00 (r, x_|_) contours because the velocity of the bubble walls are 
smaller than in the ideal e; cLS CcLSG while the ejected cooler hadronic matter tends to have 
higher flow velocity (see p4| ). 



■^Bubble formation is a natural characteristic of a mixed phase, the above foam, however, is induced by 
the initial inhomogeneities and the subsequent collective motion. 
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In Figure 14, a cut along the x— axis provides a close-up of the complex nature of the 
sharp structures of mixed phase matter that emerge as well as of the chaotic transverse 
velocity fields in between them. It also shows the energy density scales corresponding to 
the shaded contours in Figs. 12 and 13. In this view, the sharp initial inhomogeneities due 
to mini-jets superimposed on the smooth beam-jet background and the initial turbulent 
velocity field are particularly clearly revealed. 

The evolution of the transverse energy in this event is shown in Fig. 15. Unlike in 
the static examples discussed before, the initial dE^/dyd(p is not azimuthally symmetric in 
this case because of the initial turbulent velocity field. Of course, the azimuthal angles of 
the bumps and valleys vary from event to event. That is why correlation functions must 
be studied experimentally. However, the evolution of dE±/dyd<p in this event reveals the 
general tendency of inhomogeneous initial conditions to evolve into multiple azimuthally 
directed flow structures. Comparing the ideal gas and first order transition cases, we see 
again that the former leads to a lower average final transverse energy than the latter due to 
extra work done by the ideal gas upon longitudinal expansion. 

In Figure 16 we show dE±/dyd(f> averaged over 50 events and Cet(A0) for such turbu- 
lent initial conditions. The event-averaged dEj_/dyd(fi remains approximately azimuthally 
symmetric as required. Initially, there is a small (< 1%) azimuthal auto-correlation that is 
induced when the HIJING parton data are coarse-grained into 1 fm 3 fluid cells. The initial 
state correlation also includes the small fluctuating dipole contribution arising from the fact 
that P±a 7^ for a finite number of mini jets in the central rapidity slice. At later times, 
however, a 3% auto-correlation develops in the ideal fluid case and approximately 2% in the 
first order transition case. The magnitude of the correlations for both equations of state is 
evidently small. 

In Figure 17, the dependence of the induced E± correlations on the soft background as 
well as on the mini-jet scale parameter p is shown. The HIJING default case from Fig. 
16 corresponds to the solid curves. In the ideal gas case (left panels), the auto-correlation 
reaches 6% if the background is reduced by a factor of two (dashed curve) to dE^/dy = 250 
GeV as in Figs. 2-4. On the other hand, with po = 2 GeV/c fixed but the soft background 
increased by a factor of two (dotted curve) to dE s ° /dy = 1 TeV, the final auto-correlation 
is reduced by around a factor of two relative to the default case. Finally, if p$ is decreased to 
1 GeV/c but dE s °^ 1 /dy = 0.5 TeV is kept fixed (dash-dotted curve), the correlation hardly 
changes relative to the corresponding p = 2 GeV/c case (solid). The wiggle in this last 
case is due to the more limited statistics (20 events) available for this average. We conclude 
that in the ideal e the collective azimuthal anisotropies are approximately linearly 

dependent on the level of the soft beam-jet component. 

The initial state correlation function, however, also depends on the po and soft background 
scales. In the lower panels, the ratio of the final correlation function to the initial one is 
shown in the restricted A0 < 7r/2 interval that avoids the artificial pole created at the point 
where the initial correlation function crosses zero. It is interesting to note that this ratio 
for the three po = 2 GeV/c curves is practically independent of the background level. On 
the other hand, the p — 1 GeV/c ratio is significantly larger. Thus, while the absolute 
magnitude of the correlation function depends roughly linearly on the soft background level, 
the dynamical enhancement of the initial state correlations in the ideal gas case peaks near 
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4-5 at A0 = approximately independent of that background. 

While the absolute collective signature as measured by the ^-correlation function is 
small, it is quite significant compared to the initial state correlations (Fig. 16) and the 
numerical accuracy (Fig. 7). The transverse energy correlation function is of course only one 
way to search for azimuthally asymmetric collective flow phenomena. The power spectrum, 
E±(m) = (J d<p e~ mr ^ E±(4>)), wavelet analysis, and factorial moment fluctuation analysis 
may provide more sensitive probes of the induced collectivity that is apparent in the ratio 
curves in Fig. 17. 

In the case of the first order phase transition, the Et correlations are significantly sup- 
pressed relative to the ideal g £IS CclS6 and appear to be much more insensitive to the back- 
ground level. Comparing the solid and dotted curves indicates, however, a stronger depen- 
dence on the p mini-jet scale. For p = 1 GeV/c and dE s ± /dy = 1 TeV, the azimuthal 
anisotropies fall below 1%. On the other hand, in the ratio of final to initial correlation 
functions, the largest enhancement occurs for the po = 1 GeV/c case. The ratios also show 
a qualitative shoulder feature in the first order transition case for which we have not found 
a simple explanation. The level of collectivity is, however, probably too small to see such 
structures. 

The suppression of collective flow phenomena in the case of a first order transition is due 
to the vanishing of the speed of sound over a large interval of energy densities, e# < e < eg. 
This is also manifest in the smaller reduction of the initial dE±/dy due to longitudinal 
expansion relative to the expansion with an ideal gas equation of state. The above results 
can be understood as a consequence of a rather general feature of evolution with any equation 
of state that possesses a dip of c 2 s = dp/de in a finite range of energy densities. In the case of 
a first order transition, c 2 s = in the mixed phase, and pressure gradients driving collective 
flow phenomena are strongly suppressed. As emphasized in ||, even a continuous cross-over 
transition may feature such a minimum if the cross-over temperature region is not too broad. 
For realistic AT/T C ~ 0.1, the softening of the equation of state is sufficiently strong that 
the time-delay signature discussed in the next section should still be observable. The same 
physics of softening is expected to lead to a suppression of directed transverse flow phenomena 



at much lower AGS energies [f|T] and also to the suppression |42| of d(p±) / d(dN w / dy) in the 
mixed phase region. In the turbulent-glue scenario, the suppression of pressure gradients in 
the cross-over region of the equation of state has the observable consequence of reducing the 
azimuthally asymmetric collective phenomena. 

It is indeed curious that many "barometric" signatures of QGP formation involve the 
suppression of collective behaviour that would otherwise naturally arise in ideal gas systems. 
This makes the search for signatures of the QGP more difficult because ordinary dissipative 
effects due to viscocity, thermal and colour conductivity etc. work in the same direction. It 
is only through the careful systematic studies of all observables as a function of beam energy 
and nuclear size that one can hope to unravel interesting threshold-type behaviour caused 
by the passage through a mixed-phase region from the suppression of collective phenomena 
due to less interesting dissipative dynamics. 
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4 Robustness of the Time-Delay Signature 



Since it is the reduction of collective observables in the evolution of a QGP that signals 
rapid cross-over regions in the equation of state, it is important to find as many correlated 
signatures as possible in the search for that phase of matter. As repeatedly stressed above, 
one of the generic consequences of hydrodynamics with an equation of state that has a soft 



that signature of stalled dynamics is fortunately robust to an increase in the width of the 
cross-over region. In this section, we want to demonstrate in more detail the robust character 
of that observable even to the much more unfavourable turbulent initial conditions discussed 
above. 

In Figure 18, the evolution of the mean energy density is shown, averaged over the inner 
r± < 3 fm core of the plasma. The average is again over 50 events for each equation of 
state. At r = T t h = 0.5 fm/c, the mean central energy density is approximately 16 GeV/fm 3 
for the turbulent ensemble with p Q = 2 GeV/c and dE± /dy = 0.5 TeV. The solid curve, 
for the case of a first order transition and 3+1-dimensional expansion, should be compared 
to the thick dashed curve for the case of an ideal gas equation of state. In addition, the 
light dashed and dash-dotted curves are shown for comparison. They correspond to ideal 
one-dimensional Bjorken expansion with transverse expansion neglected. The r _1 curve rep- 
resents pure longitudinal expansion without work, p = 0, while the r _4//3 curve corresponds 
to the ideal gas case. The 3+1-dimensional ideal gas evolution starts to deviate from the 
one-dimensional case after a short time ~ 2 fm/c due to rapid radial expansion. In the first 
order transition case, the mean energy density follows the ideal one-dimensional Bjorken 
curve up to ~ 6 fm/c because the transverse expansion is stalled. The freeze-out occurs near 
e ~ tn — 0.1 GeV/fm 3 . It is clear from Fig. 18 that the freeze-out time is approximately 
twice as long in the case of a phase transition to the QGP. 

Especially important for the R ou t/Rside signature || 0, ||] of the QGP is that the transverse 
coordinate distribution in the first order transition case remains more compact even though 
it takes a longer time for the system to freeze-out. This can be seen clearly in Figs. 12 and 
13 which, together with Fig. 18, confirm that the space-time geometry of freeze-out remains 
so different in the two cases. 

Figure 19 emphasizes the strongly inhomogeneous character of expansion in the turbulent- 
glue scenario. The thick solid and dashed curves are the same as in Fig. 18. The two thin 
curves show the magnitude of the rms fluctuations of the energy density in the central 
core for the first order transition case. The initial fluctuations are already large. However, 
those fluctuations grow rapidly as the system passes through the mixed phase and the foam 
structure in Fig. 13 develops. This shows that freeze-out is not a fixed-time event, but is 
spread over the entire evolution. The most essential aspect for the time-delay signal is that 
the freeze-out time duration is large relative to the transverse coordinate dispersion || [7[ H . 

The actual computation of the pion interference pattern from the turbulent evolution is 
beyond the scope of the present study. A substantial generalization of the already computa- 
tionally demanding methods used in Ref. |5| will be required. However, based on that work, 
we expect the general interference pattern to reflect well the underlying, time-delayed freeze- 
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out geometry. Therefore, the time-delay signature survives even if the initial conditions are 
as inhomogeneous as in the turbulent-glue scenario. 

The detailed interference pattern can be expected, on the other hand, to differ from the 
patterns produced from the evolution of homogeneous initial conditions due to the overlap- 
ping bubble wall geometries found in Fig. 13. The Fourier transform of such hadronic foam is 
bound to contain extra structure since multiple distance scales (shell thickness and radii, and 
relative separations) enter the coordinate space distribution. It would be interesting in the 
future to explore such novel interference patterns that would be specific to inhomogeneous 
geometries. 



5 Enhanced Radiance of Hot Spots 

A specific consequence of hot-spot formation is enhanced radiance of hard probes. For 
instance, the (invariant) rate for emitting photons with momentum k from a fluid element 
consisting of quarks and gluons with temperature T and 4-velocity w M is given by ||3| 



It was moreover shown in [|3| that this rate is approximately the same if the fluid element 



consists of hadrons instead of QGP. To estimate photon radiation in the turbulent-glue 
scenario as compared to that from a homogeneous Bjorken cylinder, we integrate numerically 



the rate (18) over the photon emission angle, over (proper) time and transverse coordinates, 
and evaluate it at y = rj = (rj = \n[(t + z)/(t — z)]/2 being the space-time rapidity). 

The final differential photon spectrum diV 7 /dT]dyk±dk_L\ y=ri= o is shown in Fig. 20 av- 
eraged over 10 IIIJING events (thick lines) and for a homogeneous Bjorken cylinder (thin 
lines) of radius R = 4 fm with the same initial (average) energy as the HIJING events. We 
note that even a single event produces a quite similar spectrum as the 10-event average. As 
one expects, the higher temperatures of hot spots in the turbulent-glue scenario lead to an 
enhancement of the exponential tail of the photon spectrum. 

Comparing the yield in the case of a phase transition (dashed lines) to that in the case of 
an ideal gas equation of state (solid lines), one observes that the longer lifetime of the system 
in the case of a phase transition leads to enhanced radiation of photons with small as well 
as with large momenta in the case of the turbulent-glue scenario, while for the homogeneous 
Bjorken cylinder, only radiation of photons with small momenta is significantly enhanced. 
The higher fluid velocities in the case of the first order transition scenario as well as reheating 
in the shocks created in the expansion of hot spots seem to be responsible for the stronger 
population of the high-momentum tail of the spectrum. Correlations of direct gammas with 
azimuthally directed transverse shocks should also be looked for. Other probes such as 
strangeness enhancement may also be correlated with transverse shocks. 
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6 Summary 



In this paper we studied several possible consequences of initial state inhomogeneities and 
turbulence in quark-gluon plasmas arising from multiple mini-jet production in ultrarel- 
ativistic nuclear collisions. The HIJING model was used to generate the initial mini-jet 
configurations and to provide an estimate of the soft beam-jet background. The ensemble 
of initial conditions was found to exhibit a wide spectrum of fluctuations not only of the 
initial energy density distribution but also of the initial transverse velocity field. The fluc- 
tuations are large if the equilibration time of the plasma is short, ~ 0.5 fm/c, as suggested 
by recent estimates of radiation energy loss. We refer to this type of initial conditions as the 
turbulent-glue scenario to contrast it with the more conventional hot-glue scenario which 
assumes cylindrical symmetry, homogeneity, and quiescence of the initial plasma. 

We assumed the validity of non-dissipative hydrodynamics in order to assess the most 
optimistic observable consequences of such unfavourable initial conditions. We studied three 
observables that could serve as diagnostic tools of such plasmas. 

First, we showed that new types of azimuthally asymmetric collective flow patterns (vol- 
canoes) could arise in inhomogeneous geometries. At the event-by-event level they could be 
observed by looking for spikes in the transverse energy distribution, dE±/dyd<f). However, 
the collective effects are difficult to identify because of large uncorrected fluctuations asso- 
ciated with the turbulent nature of the initial conditions. Those uncorrelated fluctuations 
can be averaged out by studying instead the azimuthal correlation function of the transverse 
energy. The remaining dynamical correlations are small but sensitive to the plasma equation 
of state. In general, evolution with an equation of state featuring a minimum of the speed of 
sound in the energy density range of interest reduces all collective flow phenomena relative 
to the ideal e; cLS CcLSG. For example, the overall reduction of the initial transverse energy due 
to work associated with longitudinal expansion is maximal for the ideal g as case. 

Second, we discussed the time-delay signature of a QGP transition. We found that, as 
in the conventional hot-glue scenario, the evolution is stalled if there is a soft region of the 
equation of state. It is remarkable that this phenomenon survives not only if the QGP 
equation of state only features a smooth cross-over instead of a first order transition [|], Q , 
but also if the ensemble of initial conditions is as complex as in the turbulent glue-scenario. 
Meson interferometry || |7|, || appears therefore to be one of the most robust diagnostic tools 
in the search for the QGP transition. 

Finally, we considered briefly the effect of hot spots on hard probes. The thermal con- 
tribution to probes such as direct photons is obviously exponentially sensitive to the local 
temperature. We showed that the turbulent-glue scenario can enhance by more than an 
order of magnitude the high-Zc^ tails for both equations of state. Other hard probes and 
even softer ones such as strangeness production can be expected to be correlated and en- 
hanced due to hot spots and collective transverse shock phenomena. However, to identify 
any enhanced yields with hot-spot formation will require subtraction of other pre-equilibrium 
contributions to those yields. 

The above scenario represents in many ways the most optimistic point of view for sig- 
natures of QGP formation. We neglected all dissipative effects that tend to dampen col- 
lective behaviour of the system and decrease the freeze-out time. The thin expansion shell 
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structures found may be diffused considerably by such effects. Chemical equilibrium may 
not be maintained through the transition point. In general, chemical rate equations would 
have to supplement the hydrodynamic equations. Also, we assumed Bjorken boundary con- 
ditions and hence neglected fluctuations along the longitudinal direction. In the general 
3+1-dimensional case, fluctuations will produce hot-spot droplets instead of the cylindrical 
structures studied here. That extra dimension decreases the overlap phase space of expand- 
ing inhomogeneities and therefore reduces the pattern of collective transverse shocks. We 
have also assumed a most conservative homogeneous beam-jet background underneath the 
turbulent mini-jet plasma. In confronting data eventually, these and other complications 
will have to be considered in more detail in future studies. 
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Appendix: The Transverse Energy Flow 

The local thermal phase space density for a mixture of massless fermions and bosons is 
fd d e N _ 1 f g B g F \_ fp-u\ 

where T(x) and u^(x) are the local temperature and fluid 4-velocity, and Qb^f) are the 
number of active massless Boson (Fermion) helicity states. For example, in the mixed phase 
qb = 16r<g + 3(1 — tq) if the fraction of matter in the plasma phase is tq. The fluid 
energy-momentum tensor is related to / via 

= 75^* Or 1 ) • < 20 > 

For longitudinal boost-invariant (Bjorken) boundary conditions, T = T(r, Xj_), depends only 
on the transverse coordinates and the proper time r = y/t 2 — z 2 . The longitudinal fluid 
velocity is determined by the space-time rapidity variable rj = ln[(t + z)/(t — z)]/2, via 
v z = u z /u° = tanhr?. If the transverse flow velocity at z = (rj = v z = 0) is specified as 
vj_ = tanhy^(cos0^, sin0j_) in terms of the local transverse flow rapidity y±(r, xj_) and its 
azimuthal direction by the angle 0j_(r, xj_), then the fluid 4-velocity is given by 

u^{r, xj_, rj) = (cosh rj cosh y±, sinhy_|_ cos (f)±, smhy± sin cj)±, smlir? cosh y±) 

= 7_i_(cosh77, v± cos0^, v± sin sinh?]) , (21) 
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where 7_i_ = cosh?/j_ and v± = tanhy±. Note that if the transverse flow velocity vanishes, 
y± = 0, then we recover the familiar Bjorken flow velocity field u" = x^/t. 

It is convenient to express p" in terms of longitudinal rapidity and cylindrical transverse 
momentum variables, such that p" = p±(coshy, cos0,sin0,sinh?/). Therefore, the argument 
of the thermal distribution h above reduces to 

p ■ u/T = p_i_7_i_(cosh(y — 77) — v j_ cos(0 — 0j_))/T (22) 

It is linearly dependent on p± because we assume that all components of the fluid are 
massless. 

The energy density at z — 0, i.e. rj = 0, is then given by £ = T 00 via 

£(r,x ± ,0)= — / d?/cosh 2 W — V - — / dxx 3 h(x). (23) 

\1± J Jo (cosh y — v± cos <p) Jo 

Note that 

dxx 3 h(x) = —{g B + 7g F /S)— . (24) 

47T 60 

When this last factor is combined with T 4 we obtain the local proper energy density e(r, x^) 
/(47r). In the limit where v±(t, Xj_) = 0, the remaining two integrals then yield a factor Air 
and thus eq. (^3|) reduces to the expected 

lim £ (r.xj.,0) = (<? B + 7^/8) |- T 4 (r,x ± ) = e(r,x ± ) . (25) 

When the local transverse flow velocity does not vanish, the last two integrals are still 
analytic and implement the correct Lorentz transformation of the energy density in the ideal 
fluid limit with 

£(r,x ± ,0) = i(4 7 i-l)e(r,x ± ) . (26) 

(Recall that for an ideal fluid T" v = (e + p)u»u v - pg^ v .) 

If the fluid is frozen out at a fixed proper time, t/, then the invariant total particle 
distribution is obtained by integrating over all fluid cells via the Cooper-Frye formula 



dN dN 

P 



d 3 p dy d 2 pj 



/ d^Pn /(t/,x,p) 



T f p± 



J d 2 x ± J d V coshfo -y)h (^) , (27) 



where we used the 3-volume element 

da" = ^T^^i^l^dV = (coshr/,0,0,sinhr/) r f dr/d 2 x ± (28) 
or] ox± oy± 

given x M (r, xj_, r?) = (r coshr?, x_|_, rsinhr/). Note that dA r /dyd 2 pj_ is of course independent 
of t/ given our boundary conditions, but if ?/i(r,xi) ^ it could be highly azimuthally 
asymmetric. 
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Finally, the transverse energy distribution at proper time 17 and y = is given by 



dE ± 
dyd(p 



dN 



dp± p ± 



dyd 2 p_L 

47T7j 



Tj / d 2 x_|_ / dr/ cosh 77 / dp± p\h 



p ■ u 



di] 



cosh 77 



(cosh 77 — v± cos(0 — 4>±)) 4 



J d 2 x ± r / £:(r / ,x_ L ) A(<f> - 0_l(t/, xj, v x (r f , x ± )) 



where the smeared-out azimuthal angle delta function is given by 



A(<M 



di] 



cosh?] 



47T7 4 (47 2 — 1) J -00 (cosh?7 — v cos0) 4 

1 (1 -v 2 ) 3 jw(2w 2 + 13) 6 + 24w 2 
47 3 + t; 2 I (l-«; 2 )3 + (1- W 2 )V2 arctan 



'1 +w 
1 — w 



i 



(29) 



(30) 



where w = vcos<p. In the v = limit, A(0, 0) = 1/8 which is it/A times smaller than the 
naive l/(27r) uniform distribution expected in the non-thermal Bjorken limit. In the v — > 1 
limit, we can use the integral H 



2tt 



(a — &cos0) 4 
to show that J d0A(0, t> — ► 1) = 1, i.e., 



7T 



q(2a 2 + 3fe 2 ) 

( a 2 _ {,2)7/2 



(31) 



lim A(0, 

V— >1 



f32) 
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Figure Captions: 



Fig. 1: HIJING1.3 model || predictions for central (6 = 0) Au + Au collisions at y/s = 200 
AGeV are shown for the mini-jet cut-off scale p = 1 GeV/c (solid) and 2 GeV/c (dashed). 
The rapidity density of mini-jet gluons including initial and final state radiation is shown 
in part (a). Shadowing and jet quenching are not included in these calculations. The gluon 
transverse energy per unit rapidity is shown in part (b). In part (c) the initial transverse 
momentum distribution of the mini-jet gluons at y — is shown. The final hadronic trans- 
verse energy distribution including both mini-jet and beam-jet fragmentation is shown in (d). 

Fig. 2: (a) Hot spots in Au + Au (\/s = 200 AGeV, b = 0) at the thermalization time 
t = T t h = 0.5 fm/c are seen as spikes in the local energy density, £(7^, x_i_), as a function of 
the transverse coordinate, xj_, in the z = plane. The transverse and longitudinal resolution 
scales are fixed to be Ar_|_ = 1 fm and Ay = 1. The mini-jet scale here is po = 2 GeV/c. 
(b) The chaotic initial transverse momentum field, Mj_(r t ^, Xj_), is represented here by ar- 
rows, (c) The spectrum of proper energy density fluctuations at in the central region 
with r± < 4 fm is shown. The histogram is an average of 200 HIJING events and includes 
a soft gluon component due to beam-jet fragmentation, (d) The spectrum of effective gluon 
temperatures in the central region corresponding to (c) is shown. 

Fig. 3 : A surface plot of the local energy density £(r^,xj_), averaged over 200 events, 
shows the smooth, azimuthally symmetric form usually assumed in the hot-glue scenario. 
The shaded contour on the top part is a 2-dimensional projection of that surface plot. In 
parts (b)-(d), three separate HIJING events show that the fluctuations in Fig. 2a are typical. 
The shaded contours illustrate the complex nature of the azimuthally asymmetric inhomo- 
geneities produced by the mini-jet mechanism. 

Fig. 4 : The gluon number densities corresponding to the same events as in Fig. 3 show that 
the hot spots in Fig. 3 are not caused by isolated jets but occur as accidental coincidence of 
several gluons in the same region. 

Fig. 5: Evolution of a cylindrically symmetric Gaussian hotspot of radius 2 fm and peak 
energy density 8 = 30 GeV/fm 3 . Shown are (calculational frame) energy density profiles for 
(a) the initial situation at time to — 0.5 fm/c, (b) at time t = 14.9 fm/c after evolving with 
the standard antidiffusion fluxes and with the ideal gas equation of state, (d) at the same 
time after evolving with antidiffusion fluxes reduced to 70% of the standard value, and (c) 
at the same time after an evolution with the equation of state with a phase transition and 
reduced antidiffusion. (b)-(d) show the profiles only for positive x, y, the other quadrants 
are symmetric. 

Fig. 6: (a) Energy density profiles along the x— axis for times t = r + 1.6 n fm/c, r = 0.5 
fm/c, n = 0,1,..., 9, calculated with the ideal gas equation of state. Full lines are results 
from the RHLLE, dots are cell values for the SHASTA run. (b) Energy density profiles for 
times t = tq + 3.2n fm/c, n = 0, 1, ..,4 and t = 14.9 fm/c, calculated with the equation of 
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state with a phase transition. 



Fig. 7 : The evolution of numerical errors due to the cartesian grid is shown for the expan- 
sion of a single Gaussian initial condition. The transverse energy fluctuations are initially 
azimuthally symmetric but develop 10 % systematic asymmetries in the case of a first order 
transition (upper right panel). The maximal E± correlation induced by the grid is, however, 
only on the order of 0.1% (lower right panel). In the ideal gas case, errors are about an order 
of magnitude smaller (left panels) at the same time. 

Fig. 8 : The evolution of the energy density contours from the expansion of two adja- 
cent Gaussian hot spots is shown for the case of an ideal gas equation of state. Note the 
development of two back-to-back transverse jets from the shock zone. The time step between 
different frames is 0.8 fm/c, the initial time is 0.5 fm/c. 

Fig. 9 : Azimuthally asymmetric transverse shocks produced during the evolution of the 
inhomogeneous initial condition in Fig. 8 are evident as spikes in the transverse energy dis- 
tribution (upper panels) as well as enhanced Ej_ correlations (lower panels). The initial 
dE±/dyd(f) and Cet are straight lines, due to the absence of collective motion in the initial 
conditions. Note that a first order transition (right panels) inhibits the formation of such 
transverse shocks because the expansion is stalled in that case (see Fig. 10). In both cases, 
the dynamically produced asymmetries are much larger than the numerical errors from Fig. 7. 

Fig. 10: The evolution of the same initial condition as in Fig. 8 but assuming here the 
Bag model equation of state. The stalled expansion due to the vanishing of the speed of 
sound is revealed using larger time steps (1.6 fm/c) between frames relative to Fig. 8. Note 
that the expansion bubbles are much thinner as is the shock zone. 

Fig. 11: The evolution of the energy density £ along the a;— axis of Figs. 8 and 9 shows 
clearly the difference between shell structures formed in the ideal gas vs. the first order tran- 
sition cases. Lower graphs show the evolution of the v x (x, 0) transverse velocity field that 
vanishes initially. The time steps shown for the ideal case are r = 0.5, 2.1, 3.7, 5.3, 7.7 fm/c, 
while in the first order transition case they are r = 0.5, 3.7, 6.9, 10.1, 14.9 fm/c. 

Fig. 12: The evolution of the energy density S for a typical HIJING event assuming an 
ideal gas equation state. The (solid) contour indicates the hadronization surface, where the 
proper energy density has the value e# — 0.1 GeV/fm 3 . The time steps are indicated on 
the top right of each frame. Darker shaded regions correspond to hot spots while lighter 
shaded regions to cooler domains. See Fig. 14 below for a slice along the x— axis through 
these profiles that indicates the energy density scale. 

Fig. 13: The evolution of the same HIJING event as in Fig. 12, but with the equation 
of state with a first order transition. Note that the time steps are twice as large as in Fig. 
12 due to the stalled expansion with this equation of state. Also note the formation of thin 
bubble walls of mixed phase matter as in Fig. 10. The contour of e = e# shows that hadronic 
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foam-like structures develop due to the initial state inhomogeneities. 

Fig. 14: A slice along the x— axis through selected frames in Figs. 12 and 13 indicates the 
magnitude of the energy densities in the structures formed during evolution. The v x (x,0) 
collective flow velocity field does not vanish at T t h as in Fig. 11 because of the initial tur- 
bulence induced by mini-jet formation. The curves are plotted at the same times as in Fig. 11. 

Fig. 15: The evolution of the thermal smeared transverse energy distribution corresponding 
to the event in Figs. 12 and 13. Note that the initial azimuthal anisotropy reflects the turbu- 
lent initial condition. At freeze-out, multiple narrow spikes formed from directed transverse 
shocks (volcanoes) are superimposed on top of the broader initial state anisotropy. The 
overall magnitude of the transverse energy is reduced more in the ideal case. 

Fig. 16: The event-averaged transverse energy distribution (upper panels) is approximately 
rotation invariant within the finite statistics. However, a clear enhancement of the transverse 
energy correlations is observed (lower panels). The initial state correlations (dashed curves) 
are typically less than 1% and the final dynamical correlation due to transverse shocks (solid 
curves) are strongest in the ideal case (lower left panel). 

Fig. 17: The dependence of the E± correlation function on the soft beam-jet background 
and the mini-jet po scale is shown. The solid curves correspond to the default HIJING 
p = 2 GeV/c and E s ° ft = dE s ° ft /dy = 0.5 TeV shown in Fig. 16. For the ideal gas case 
(left panels) the dashed and dotted curves correspond to p = 2 GeV/c and E s °^ = 0.25, 1.0 
TeV, respectively. For the first order transition case (right panels) the dashed curves corre- 
spond to po = 2 GeV/c and E s °* 1 = 0.25 TeV and the dotted curves to p — 1 GeV/c and 
£soft _ y q TeV, respectively. The dash-dotted curve corresponds to po — 1 GeV/c and 
£soft _ q 5 TeV. The lower panels show the ratio of the final to initial correlation functions 
in the restricted A0 < tc/2 domain. 

Fig. 18: The evolution of the central core (r± < 3 fm) mean energy density (S) is compared 
for 3+1-dimensional (thick solid and dashed curves) and one-dimensional Bjorken expan- 
sion (thin dashed and dash-dotted curves) and different equations of state (thick solid: first 
order transition, thick dashed and thin dash-dotted: ideal gas, thin dashed: p — 0). 

Fig. 19: The evolution of the mean and rms fluctuations of the central core energy den- 
sity (£) is shown for the first order transition case. The dashed curve is the same as in Fig. 18. 

Fig. 20: Photon transverse momentum spectra at y = rj = 0. Thick lines are the aver- 
age spectra calculated from the hydrodynamical evolution of 10 (turbulent) HIJING events, 
thin lines are for the evolution of a homogeneous Bjorken cylinder of radius R = 4 fm with 
the same initial energy content. Dashed lines are obtained assuming the equation of state 
with phase transition, full lines are for the ideal gas equation of state. 
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